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Abstract. We study the transport properties of a one-dimensional hard-core bosonic 
lattice gas coupled to two particle reservoirs at different chemical potentials which 
generate a current flow through the system. In particular, the influence of random 
fluctuations of the underlying lattice on the stationary-state properties is investigated. 
We show analytically that the steady-state density presents a linear profile. The local 
steady-state current obeys the Fourier law j — ~k(t)Vti where r is a typical timescale 
of the lattice fluctuations and Vn the density gradient imposed by the reservoirs. 



1. Introduction 

The transport properties of energy or particles in small quantum systems are an 
important topic in nonequilibrium statistical dynamics. In particular, the transition 
between ballistic and diffusive transport is at the centre of many investigations, 
attempting to understand from a microscopical point of view the emergence of the 
celebrated Fourier law [T] [2] . With the development of nanoscale technologies, it is now 
becoming possible to design proper experiments that can potentially test theoretical 
predictions for small quantum systems. The most promising possibilities certainly come 
from the optical lattice community. For example, optical lattices can now be used to 
experimentally generate one-dimensional (ID) bosonic systems [31 HIE] which have been 
studied theoretically for many years [6j[7l|8]. In the large scattering-length limit and at 
low densities, ultracold bosons effectively behave as impenetrable particles [9], namely 
as hard-core bosons, thus realising the Tonks-Girardeau model [SI EJ. Experiments 
on such ID hard-core bosons have been performed with Rubidium atoms within both 
continuum [5] and lattice contexts [10J . For a recent review of developments in ultracold 
gases and optical lattices, see [TT] . 

In this letter we present the transport properties of a ID-lattice hard-core bosonic 
gas driven out of equilibrium by the interaction at its boundaries with two external 
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reservoirs that induce a particle-current flow through the system. We remark here that 
similar studies have been performed in [12j DSl EH] • The particular focus of the present 
work is the influence of lattice fluctuations, which may either be induced artificially or 
be inherent to an experimental set-up, on the transport properties 

2. Model 

The Hamiltonian associated with the hard-core boson model on a linear optical lattice 
of N sites is given by 

N N-1 

Hs = J2 hl + J2 Vl - 

i=i i=i 

Here the on-site one-particle Hamiltonian is 

hi = eb+k = em , (2.2) 

with a site-independent chemical potential e coupled to the local occupation number 
ni = b^bi, while the hopping potential is 

V t = -h [b+b l+1 + b+ +1 k] (2.3) 

where the hopping rate t\ may depend on the position in the trap. The creation and 
annihilation operators satisfy the usual bosonic commutation relations on different sites, 
[&,,&+] = [b h b v ] = [bf,b$] = 0, while the hard -core constraint is implemented by the 
additional conditions bf = (bf) 2 = and {bi, b^} = 1 preventing more than single 
occupancy of sites. Notice that through the transformation b + = (o~ x + io~ y )/2, where 
a x,v are the usual Pauli matrices, the hard-core boson Hamiltonian is exactly mapped 
onto the XX quantum spin chain which, in recent years, has been studied extensively 
in a nonequilibrium context [T5] . 

The bosonic gas inside the trap is coupled at its left and right boundaries to 
ideal (non-interacting) hard-core bosonic reservoirs set at different densities and 
ur, described by the single-particle density matrices 

p L ,R=\l}n L , R (M + \0)(l-n LtR )(0\ (2.4) 

where the labels L and R stand for the left and right reservoirs respectively and |0), |1) 
are the associated vacuum and one-particle states. The interaction with the reservoirs 
is implemented via a discrete-time repeated interaction scheme which in the continuum 
limit leads to a Markovian Lindblad dynamics |16j . See [17] for a useful discussion on 
the possible failure of Lindblad dynamics in the description of stationary nonequilibrium 
properties and the importance of neglecting the internal couplings. Within the discrete 
process, at a given time t only one left reservoir particle and one right reservoir particle, 
in state pi and p R respectively, interact with the system. These particles interact for 
a time r through the hopping potential V and V^. After the interaction, i.e., at time 
t + r, the system state ps = Tr^{p}, obtained after tracing out the environment degrees 
of freedom corresponding to the left and right reservoirs, is given by 

p s (t + t) = Tr L<R {Uj {p L ® p s {t) ® PR ) U^} . (2.5) 
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Figure 1. Sketch of the time-evolution of the fluctuating optical lattice. Here the left 
reservoir is empty while the right one has full occupancy. Local fluctuations, underlined 
by a short straight line, enhance the hopping rate locally. 



Here Uj = e 1tHt with the total relevant Hamiltonian given by = H$ + Vq + Vjv + 
ho + h^+i where /io,tv+i are the one-particle Hamiltonians of the reservoirs. Notice here 



that Ht is of the form (2.1) with N + 2 sites. The process is then repeated with new 



reservoir particles such that (2.5) is iterated further. The net effect of the process, for 



every timestep r, is that a boson can be injected into the trap or escape from it. For 
example, the extreme limit = and = 1 describes the injection of bosons from 
the right of the trap and their escape to the left, i.e., in this case escape to the right 
and injection from the left are forbidden. 

As mentioned in the introduction, we consider the effect of fluctuations of the optical 
lattice which may be induced by some underlying physical process, e.g., vibrations of 
mirrors, presence of impurities. We simply model this randomness by allowing that 
the hopping rates t\ VZ = 0, ...,N + 1 fluctuate in time within a typical timescale 
t f. In the following we assume Tf ~ r. During the time-evolution each hopping rate 
follows a stochastic trajectory, see figure [TJ which is governed by some known probability 
distribution. 



3. Dynamics 

Given an initial equilibrium system state, ps(0), we start the dynamics by iterating 



(2.5) with the evolution operator Ui following the fluctuations of the hopping rates. 



Instead of solving directly the dynamical equation for the density matrix, we study the 
time-evolution of correlation functions. Moreover, due to the free-fermionic structure of 
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the model after a Jordan- Wigner transformation [18 



r« = 4 = + „ = o,...,* + i (3.i) 

r K+2+ , = -iB, = -ie"Ej-.«/(6, - 6+) 

where the Ts are Majorana real (Clifford) operators satisfying = F and {r^Tj} = 
25ij, thanks to Wick's theorem, one can express all physical observables in terms of the 
two-point correlation functions 

[G(% fe = ^Tr{[r fe ,r^(t)} • (3-2) 

In the Heisenberg picture, the time evolution of the Majorana field T, generated by Ht, 
is simply given by T(t) = e~ lt3 T(0) = R(t)T(0), where T is defined by the Hamiltonian 
in terms of the field T: = (l/4)r'TT. The matrix elements of the rotation matrix R 
are simply expressed in terms of the spectral properties of Ht, see [19] for the explicit 
forms. 

We order the = (T* E , T* s ) such that the first part, T E , is associated with the 
interacting part of the environment and the second part, Ts, with the components of 



the system. Projecting (3.2) onto the system part, one arrives at the fundamental 
dynamical equation for the system correlation matrix Gs- 

G s (t + r) = R s G s (t)R s ^ + R SE G e Rse ] ■ (3.3) 

The 2N x 2N matrix R$ is that part of the full rotation matrix R = e~ lrT = 

Re Res 
Rse Rs 

given by the lower off-diagonal block of R expressed in the basis (r^,r^). For non- 
interacting dynamics, i.e., a closed system, Rse = and the rotation matrix splits into 
a block-diagonal form where Rs and R E are the rotation matrices of the system and 



which acts on the system. The 2N x 4 rectangular matrix Rse is 



environment part respectively. In the dynamical equation (3.3), the bath properties 
enter only through the initial environment particle states, encoded in the two-point 
correlation matrix Ge- The correlation matrix Ge stays constant in time since at each 
step of the repeated interaction procedure the bath particles are replaced by fresh ones. 

4. Steady-state current and density 



In the following, we concentrate mainly on the asymptotic properties of (3.3). The 
steady state is reached exponentially with a relaxation time depending on the system 
size N. For the non-disordered situation the timescale needed to reach the steady state 
behaves as N 3 [20] • In particular, we focus our attention on the transport properties of 
the bosonic gas through the optical trap. For that we compute the density profile and 
the particle current along the chain. Since the hopping dynamics conserves particles, 
one may naturally define the particle current through the Heisenberg equation of motion 
for the density: 

h^iiHs^ij^Ji-x-Ji (4.1) 
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where Ji denotes the particle-current operator associated with the Zth bond and given 
by 

^ = = it ^6,6+ ! - . (4.2) 

The density h\ and the current J/, are easily expressed in terms of the Majorana field 
T and their expectation values ni = (hi) = Tr{hip} and ji = (Ji) = Tr{ J L p} are given 
by on-site and two-site two-point correlation functions Gs'- ni = (1 — (Gs)i,i+n)/2 and 
ji — (Gs)i,i+i- In the following the stared quantities n* and j* have to be understood 
as expectation in the steady-state. 



4-1. Non-fluctuating lattice 

If the lattice is free of any disorder, the hopping rates are uniform, i.e., t[ = t VZ. In 
this case, an exact solution of the steady state has been given in [21]. It is found that 
the density profile is flat except at sites 1 and N which are directly in contact with the 
reservoirs. The density in the flat region is given by the mean value set by the reservoirs 
n* = n — (ni+nR) /2 V/ ^ 1, N while the boundary values are n\ = n— A(t )(nR— nz,)/2 
and n* N = h + A(t )(n R — n L )/2 with a shift from the mean value h depending on the 

2 

density difference tir — ul and where A(t D ) = with 7 = t a /2. The steady-state 
current j* takes a constant value j* = —a(nR — Ul) independent of the system-size 
where a is a non-monotonous function of the hopping rate t Q , with a maximum current 
state at t a = 2. The size-independence of j* signals the ballistic nature of the transport, 
which is ultimately related to the integrability of the model, that is to the equations of 
motion of the free quasiparticles describing the system. In this case there is no finite 
conductivity k and the system obviously does not obey Fourier's law. One may notice 
that this behaviour is very similar to the behaviour observed in the classical Reider- 
Lebowitz-Lieb model of a homogeneous harmonic chain |22j in contact with stochastic 
heat baths at the boundaries. 



4-2. Fluctuating lattice 

Next we study the effect of fluctuations of the lattice parameters on the steady-state 
properties. We consider local fluctuations in the sense that within the timescale r of the 
fluctuation only one bond is affected leading to an enhancement of the local hopping 
rate from its unperturbed value t a to a larger value tp, which we choose to be 1/2. 
Moreover, we take the limit of a strongly localized lattice gas, with t„ <C 1. In other 
words, at each timestep r of the dynamics, a single bond is activated at random and 
locally the particles are exchanged with a rate to = 1/2, either within the system if 
the selected bond is a system one or with the reservoirs if the fluctuations act close to 
the boundaries. These particle exchanges are reminiscent of the well-studied symmetric 



exclusion process [23]; the dynamics (3.3) leads to a non-trivial dependence of the bulk 



density gradient on the interaction time r, as we shall now see. 
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In the weak-hopping limit t — > 0, the dynamics simplify considerably since at a 
given time step only one bond is activated. The time evolution of the system is most 
simply expressed in terms of Dirac fermions Ck = (A/~ + Bk)/2 and c\ = (Aj. — Bk)/2. 
One has in matrix form c(t + r) = e c(t) where K is the coupling matrix defining 
the Hamiltonian H T = —c + Kc in terms of the fermi field c + = (cq , cf, . . . , c^, c^ +1 ). 
Due to the fact that at a given time step only one bond is non-vanishing, suppose bond 
I conecting the sites I and 1 + 1, one has the trivial dynamics c^it + r) = e tte Ck(t) 
Wk 7^ (1,1 + 1) and ci(t + r) = e tt£ [cos |cj(i) + i sin |cj + i(t)] and c/+i(t + r) = 
e lt£ [i sin |q(£) + cos |q+i(£)] for the on bond operators. 

Consider the density gradient Si = (fty+i) — on bond /. This gradient is changed 
only if the activated bond is the Zth one or one of the nearest neighbours I — 1 and I + 1. 
From the dynamical equation (3.3) in terms of the fermi operators, if the hopping is 
enhanced on bond I then Si is mapped to 

S\ = Si cos r + ji sinr (4.3) 
and ji is mapped to 

j'l = —Si sin r + ji cos r . (4.4) 
Alternatively, if the activated bond is I ± 1, the updated density gradient satisfies 

S'i = Si - - (S l±1 cos r + ji ±1 sin r - S l±1 ) (4.5) 
while the updated current is given by 

Ji = Ji cos 2 +Pl ( 4 ' 6 ) 
where the Pi are proportional to correlations across two bonds. As mentioned before, 
under this dynamics the system relaxes exponentially towards a current carrying steady 
state. Nevertheless, we remark here that the periodicity of the dynamics implies a strong 
slowing-down of the relaxation to the steady state in the neighbourhood of r = n27r. 
Indeed, for even values of n the dynamical generator maps to the identity and for odd 
values it maps to a reflection dynamics which loses its relaxation properties. In the 
following analysis we avoid these special r values. 

In the steady state, the P terms vanish on average and the set of dynamical 
equations for the gradient density and current closes. At any given time, the last update 
on bond I has probability 1/3 to have resulted from activation on bond I, probability 
1/3 to have resulted from activation on bond I — 1, and probability 1/3 to have resulted 
from activation on bond 1 + 1. Consequently, the steady-state average (denoted by a 
star) gradient obeys 

8* (cos r — 1) + j* sinr = t](t) (4.7) 

where 7] (r) is a constant independent of the bond index I. Since the steady-state current 
j* = j* is constant in space it also follows from (4.7) that the gradient density is site- 
independent and we can thus omit the /-subscripts. The steady-state current satisfies 

1 2 r 

f = g [S* sin t + j* cos r] + -j* cos - (4.8) 
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Figure 2. Normalized steady-state density gradient as a function of the interaction 
time with a delta time-distribution on the left and an exponential one with mean r 
on the right. The full lines correspond to the analytical curves while the crosses are 
obtained numerically with a time average of the density gradient. 



which is equivalent to 

si n t 

f = —~ ~ t^* = ~<r)8* (4.9) 

3 — cos r — 2 cos | 

and defines the conductivity coefficient k(t). We now determine the constant i](t) by 
considering the boundary conditions. Remembering that the densities on the reservoir 
sites are fixed and that the boundary terms jo and jjv, which are initial correlations 
between the reservoirs and the system, vanish due to the repeated interaction scheme, 
one sees that an update on bond (between left reservoir and boundary site) gives 
5' = 7j<5o(l + cosr) whereas an update on bond 1 gives 5' = 5o~ \(b~i cosr+ji sinr — Si). 
The steady-state average Sq must therefore obey 

1 / 1 + cosr ^ 1 / 5* cosr + fsmT-5* \ 
S = 2 (/o— 2— ) + 2 [ S ° 2 J 

giving r)(r) = (cost — 1)5q . By symmetry, at the right-boundary we find 5^ = 5q. 
Noting that the reservoir density difference An = n^ — n^ = 25$ + (N — 1)6*, one finally 
gets for the bulk steady-state density gradient 

A^ + l+ 7 (r) v 1 

with the finite-size shift function (extrapolation length) 

7 (r) = 2-? «(r) (4.12) 

1 — cos r 

This analytical expression is compared with numerical simulation data in figure [2] 
obtained on chains of N = 30 spins and the agreement is seen to be excellent. 

So far, we have considered the somewhat unphysical situation where the 
enhancement of a local hopping rate always stays precisely for a time r. This hypothesis 
leads to the trigonometric form of the conductivity k(t) and shift function 7(1"). A 
more reasonable assumption would be to draw the duration of a local fluctuation from 
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a probability distribution /(r). In that case, one may average the dynamical equations 
for the current and the gradient density over the time-distribution f(r) which basically 
leads to replacing the trigonometric functions cos r, sin r and cos | by their expectations 
under /. For example, for an exponential distribution of interaction timescales with 
mean r D , one gets the shift function 7 = (2/t q )k with conductivity k = 3r ^ T 2+ 2 ) wn i cn 
is again in agreement with the numerical results, see figure 2. In the limit of short-time 
interaction r — >■ 0, conductivity diverges as r" 1 while the density gradient goes as t„. 
This leads to a linear vanishing of the steady-state current j* ~ r Q . 

If lattice fluctuations are generated diffusively, e.g., the updated bond follows 
a symmetric random walk, then analysing the dynamical equations along the same 
lines, we obtain again a linear profile of the particle density in the steady state: 
5* = Ap/ (N + 1 + 7(t)) where the shift function 7(7") depends on the precise definition 
of the random walk at the boundaries. 

At finite, but small, unperturbed hopping rates t L = t Q along the optical lattice, 
we observe numerically that the linear density profile survives. However, the density 
gradient is strongly attenuated by a function of the bulk hopping rate t Q , whose 
asymptotic behaviour is (Nt Q )~ l for large lattice sizes. Consequently, for large systems 
the steady-state density gradient behaves as 5* ~ 1/N 2 instead of the former 1/N 
behaviour. At the same time, the steady-state current is j* ~ 1/N which leads to a linear 
divergence of the conductivity coefficient with system size: k ~ j* /S* ~ N. This implies 
that the classical transport properties of the hard-core boson gas are induced by the 
optical lattice fluctuations only for sufficiently small Nt Q values. In the thermodynamical 
limit, N — > 00, the system shows a ballistic transport behaviour. See [13] for a similar 
discussion. 

5. Conclusion 

We have derived analytical expressions for the steady-state conductivity and density 
profile of a nonequilibrium hard-core boson model. For a perfect optical lattice, due to 
the integrability of the model, the transport properties are anomalous with an infinite 
conductivity coefficient, reflecting the ballistic nature of the excitations. On the other 
hand, when fluctuations of the underlying lattice are present and when Nt Q is sufficiently 
small, the classical Fourier law is recovered, with a linear density profile and a finite 
conductivity coefficient depending on the lattice fluctuation properties. 

Fourier's law in nature is so wide-spread that specific choices of model parameters 
should not play a decisive role in deriving it. Consequently, we do not believe that the 
dynamical fluctuations considered in the present paper have to be the generic origin of 
normal heat conduction. In general it is believed that the basic microscopic mechanism 
leading to Fourier's law is linked to the scattering of energy carriers, inducing mixing 
properties. Indeed, since the thermal conductivity, obtained from the Green-Kubo 
formula, is given by an infinite integral of the autocorrelation function of the current 
operator, this correlator has to decay quickly enough such that the integral converges 
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even in the thermodynamic limit. In the present work we have presented a possible 
way of how Fourier's law can be established in some limiting situations by introducing 
dynamical fluctuations which somehow scatters the energy carriers. The origin of such 
fluctuations at the microscopic level can be of course very diverse, as for example a 
coupling to phonon modes or external perturbations like vibrations of the mirors in an 
optical setup. The next step of this work will be the study of the current autocorrelation 
function to compare our results for the conductivity coefficient with the Green-Kubo 
expectation, see for example [21] for a study in this spirit. 
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